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Abstract — In this paper we propose a new efficient message 
passing algorithm for decoding LDPC transmitted over a channel 
with strong phase noise. The algorithm performs approximate 
bayesian inference on a factor graph representation of the 
channel and code joint posterior. The approximate inference 
is based on an improved canonical model for the messages of 
the Sum & Product Algorithm, and a method for clustering 
the messages using the directional statistics framework. The 
proposed canonical model includes treatment for phase slips 
which can limit the performance of tracking algorithms. We 
show simulation results and complexity analysis for the proposed 
algorithm demonstrating its superiority over some of the current 
state of the art algorithms. 

Index Terms — phase noise, factor graph, Tikhonov, phase slip, 
directional statistics, moment matching 

I. Introduction 

Many satellite communication systems operating today em- 
ploy low cost upconverters or downconverters which create 
phase noise. This noise can severely limit the information rate 
of the system and pose a serious challenge for the detection 
systems. Moreover, simple solutions for phase noise tracking 
such as PLL either require low phase noise or otherwise 
require many pilot symbols which reduce the effective data 
rate. 

In the last decade we have witnessed a significant amount of 
research done on joint estimation and decoding of phase noise 
and coded information. These algorithms are based on the 
factor graph representation of the joint posterior distribution. 
This framework proposed in [1 1, allows the design of efficient 
message passing algorithms which incorporate both the code 
graph and the channel graph. The use of LDPC or Turbo 
decoders, as part of iterative message passing schemes, allows 
the receiver to operate in low SNR regions while requiring less 
pilot symbols. The best known message passing algorithm for 
phase noise channels quantizes the phase noise and performs 
good approximation of the sum & product algorithm (SPA). 
This algorithm (called DP - discrete phase in this paper) 
requires large computational resources to reach high accuracy, 
rendering it not practical for some real world applications. 

In 0, an algorithm which efficiently balances the tradeoff 
between accuracy and complexity was proposed (called BARB 
in this paper). BARB uses Tikhonov distribution parameteriza- 
tions (canonical model) for all the SPA messages concerning 
a phase node. However, the approximation as defined in 0, 



is only good when the information from the LDPC decoder is 
good (high reliability). In the first iteration the approximation 
is poor, and in fact is correct only for pilot symbols. The mes- 
sages related to the received symbols which are not pilots are 
essentially zero (no information). This inability to accurately 
approximate the messages in the first iterations causes many 
errors and can create error floor. This problem is intensified 
when using either low code rate or high code rate. In the 
first case, it is since the pilots are less significant, since then- 
energy is reduced. In the second case, the poor estimation of 
the symbols far away from the pilots cannot be overcome by 
the error correcting capacity of the code. In order to overcome 
this limitation, BARB relies on the insertion of frequent pilots 
to the transmitted block causing a reduction of the information 
rate. 

In order to improve the Tikhonov approximation, in this 
paper we suggest to avoid approximating the messages related 
to the received symbols, but to approximate only phase poste- 
rior messages by a Tikhonov distribution. A major limitation 
of the resulting Tikhonov approximation is its sensitivity to 
phase slips. If the canonical model used to approximate the 
phase posterior messages is a Tikhonov distribution, then each 
message can only have one maximum point (called Phase 
Hypotheses). Therefore, all the hypotheses (mean values of 
the Tikhonov) across the code block describe a tracking plot 
of the phase noise (called Phase Trajectory). Phase slips 
occur when the estimated phase trajectory has an ambiguity, 
which can happen if different constellation symbols have 
similar likelihoods. Since the canonical model approximates 
this trajectory using a Tikhonov distribution, information is 
lost, and finally results with tracking a wrong trajectory, a 
phenomena that resembles phase slip in PLLs. We treat this 
problem by modifying the canonical model and adding a flat 
term i.e. uniform pdf to represent the possibility that the phase 
is not represented by the tracked trajectory. The scaling of this 
uniform pdf can be viewed as the probability of a phase slip. 

In this paper, an iterative message passing algorithm is 
proposed, which uses a modified canonical model for the 
approximation of messages in the SPA. This algorithm needs 
fewer pilots and can operate in a wider range of code rates 
than BARB. The algorithm does not approximate the received 
symbols phase information, but applies a Tikhonov approxi- 
mation later, on the forward and backward recursions wherein 



the phase estimation is updated. At each such recursion step, 
the posterior probability is a mixture of Tikhonov distributions. 
We apply a clustering algorithm based on the KL divergence 
in order to select which of the components of the mixture are 
going to be approximated by the Tikhonov and the rest will 
be left over as the flat hypothesis. The modified canonical 
model approach enables us to track a phase trajectory while 
maintaining a level of confidence for the tracked trajectory. 
This approach proves to be robust to phase slips and provides 
a high level of accuracy while keeping a low computational 
load. 

An approximation of a mixture of Tikhonov pdfs to a single 
Tikhonov pdf optimal in sense of KL has not existed before, 
and required a new derivation, partially presented in this paper. 

The resulting algorithm was shown in simulation to provide 
very good performance in high phase noise level and very 
close to the performance of the optimal algorithm even when 
very few pilots are present and the code rate is high. 

II. Preliminaries 

A. Directional Statistics 

Directional statistics is a branch of mathematics which 
studies random variables defined on circles and spheres. For 
example, the probability of the wind to blow at a certain 
direction. The circular mean and variance of a circular random 
variable 8, are defined in [4|: 
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We define the operator g(9) = CMVM[/(0)] (Circular Mean 
and Variance Matching), to take a circular pdf - f(9) and create 
a Tikhonov pdf g(8) with the same circular mean and variance. 
Let g(6) be a Tikhonov distribution: 



9(0) 



27r/ (£; g ) 

The following relations are obtained: 



V>3 
h{k a ) 



^Circular (/) 



= l-a 2 c 



ircular 



(./') 



(3) 

(4) 
(5) 



An alternative formulation for the Tikhonov pdf uses a 
single complex parameter Z = ke Jft 

B. Optimal Approximation of Tikhonov Mixture 

In this section we will show that the nearest Tikhonov 
distribution to a Tikhonov mixture (in a Kullback Liebler (KL) 
sense), has its circular mean and variance matched to those of 
the mixture. The Kullback Liebler (KL) divergence [3] is a 
common information theoretic measure of similarity between 
probability distributions, which is defined as: 

D(f\\g)= P /(0)logffid0 (6) 
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a Tikhonov mixture. Then the Tikhonov distribution 17(6') 
which minimizes D(f\\g) is 
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Proof: 

We will only provide an outline of the proof due to limited 

space. Let g{6) = 27rJo(fc) 
We wish to find 
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Skipping a few algebraic steps and using (0, we get 
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where /ij = arg(z^) and ki — \z%\, parameters of the i'th 
Tikhonov distribution (mode) in the mixture f(8). This is 
basically matching the circular variance of the mixture to the 
Tikhonov distribution. 

The optimal 17(6') simultaneously needs to satisfy 
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Inserting g(9) into (fTTT i and skipping several steps we get 
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Which is basically a weighted vector sum of the expectations 
of each mode in the mixture. | ■ 

III. System Model and Previous Work 

The system model and factor graph representation as in |]2] 
are considered. For the sake of clarity, we will very briefly 
review the model and iterative algorithm. We consider the 
transmission of a sequence of complex modulation symbols 
c = (co, ci, ck_i) over an AWGN channel affected by 
carrier phase noise. We assume the symbols are drawn in- 
dependency from an MPSK constellation. The discrete-time 
baseband complex equivalent channel model at the receiver is 
given by: 

r k = c k e j9k + n k k = 0, 1, K - 1. (13) 

The phase noise stochastic model is a wiener process 

6 k = fc _! + A k (14) 

where is a real, i.i.d gaussian sequence with ~ 
JV(0,al). 

The resulting Sum & Product messages are computed by 

Pf(0k)= / Pf(0k-i)pd(0k-i)PA(0k-dk-i)d9 k -i (15) 
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where M,rfc and er 2 are the constellation order, received 
base band signal and the AWGN variance respectively. The 
messages p d ,Pb,Pd,Pf and the functions p A (4). fk(ck,O k ) 
are defined in 0. 

Due to the fact that the phase symbols are continuous 
random variables, a direct implementation of these equations 
is not possible and approximations are unavoidable. In 0, a 
Tikhonov approximation is used for all the messages in the 
SPA which leads to a very simple and fast algorithm. 

In the following sections, a new algorithm will be presented, 
which approximates the SPA messages using the directional 
statistics framework. 

IV. Approximating the SPA messages 

In El, the messages Pf(0 k ), Pb(0k) and p d (0 k ) were ap- 
proximated by Tikhonov distributions and a message passing 
algorithm was derived based on the SPA recursion equations. 
In the first iteration, when there is no information from 
the LDPC decoder, this approximation provides poor results 
since pd(0 k ) is approximated as a uniform pdf. In order to 
improve accuracy, one can suggest a different approximation 
by realizing that there is no need to approximate pd{0 k )- 

Decoding the LDPC code symbols only requires P u {c k ) 
and the phase messages act behind the scenes. Therefore, in 
the computation of ( PT31 ) and J16I . Pd(0k) can be used without 
approximation. Subsequently, only < fT3T > and (\M need to be 
approximated as Tikhonov. This modification of BARB is 
shown in figure Q] for very low phase noise variance. BARB 
estimates Pf{0 k ) as pf(6 k -i) because it has no data from 
Pd(&k-i)- Our modification approximates Pf(O k ) using the 
exact Pd(0 k -i). The results show that our modification is 
better. 

However, This approximation creates a severe performance 
degradation due to phase slips. In this section we will present 
a different canonical model for the SPA messages which is 
robust to phase slips. 

A. Modified Canonical Models 

In this paper, a new approach, more robust to phase slips, 
is presented which utilizes a modified canonical model: 

Pf(0k-i) = a k -ipf(0k-x) + (1 - ak-i)^- (19) 
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Fig. 1. Comparison between approximation methods 
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And a.k-1, Pk-i are real numbers between and 1. 
The modified canonical model can be viewed as a mixture 
of two distributions: 

1) Tikhonov: l{Zl_S) is the forward phase noise hypothe- 
ses for and \zl_ x \ is the confidence level of this 
hypotheses (inverse to the estimation variance). 

2) Uniform: representing all the other possible hypotheses 
not represented by the Tikhonov above. 

In this manner we are able to more accurately estimate the 
phase posterior and be robust to phase slip events. We will 
show the derivations only for pf, but the same applies for pb 
with a proper change in indexing. 

In the previous section we have shown that the best way, 
in the KL sense, to approximate a Tikhonov mixture using a 
single Tikhonov is to match their circular mean and variance. 

Inserting the modified canonical model ( fl9l l in dl5t . results 
in a a sum of two Tikhonov mixtures, both mixtures are in 
the order of the constellation size. 
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The modified canonical model restricts (1231 to have one 
Tikhonov distribution, thus only one phase trajectory can be 
tracked. Since the tracked phase trajectory was approximated 
using pf(0 k -i), then the next step is to process only the 
mixture (|24l . It is important to understand that not all the 
components in d24i i represent the same phase trajectory and 
there may be several hypotheses which split the tracked 



trajectory. Therefore, a selection algorithm must be employed 
which selects only one trajectory from (l24l (which may be a 
clustering of several mixture components (hypotheses)). The 
selected components are clustered to a Tikhonov distribution 
while the rest of the components of ( 1231 transfer their energy 
to the flat hypotheses. 

B. Component Selection and Clustering Algorithm 

In this section we propose an algorithm which selects 
components from the Tikhonov mixture (1231 and clusters them 
to a single Tikhonov distribution. Note that D(f\\g) is the KL 
divergence between / and g. We define the input mixture to 
the algorithm as: 
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In order to combat phase slip events, the input ( |26l > can 
be one of two mixtures: First, when the previous processed 
symbol was not a pilot, then mixture d26i > is d24l >. Second, in 
case the previous symbol was a pilot then (|26T > is (1231 . It is 
assumed that if a phase slip has occurred, then the pilot will 
correct the tracked phase trajectory, so in the second case, the 
probability of the new hypotheses is 1: au = 1. This is a 
reasonable assumption in high SNR. 

Algorithm 1 Component Selection & Clustering Algorithm 

lead <— argmaxi{ \ z ^\-i } 

idx •<— lead 

for i = 1 — > N do 

if <T D then 

idx <— [idx, i) 

end if 
end for 

p f (6k) <- CMVM(a(idx)J(idx)) 
OL k <- 5k-i J2 a(idx) 

The algorithm finds the component with the highest am- 
plitude to variance ratio (amp2var) and clusters all the other 
components "close" to it (KL sense and under a specified 
threshold). The ratio amp2var describes the ratio between 
the likelihood of the hypotheses and its estimation error. The 
assumption is that all the clustered components represent the 
same phase trajectory. The likelihood &k is computed using the 
weights of all the clustered mixture components. Therefore, 
the modified canonical model can be viewed as a weighted 
sum of mixture components which better estimates the phase 
posterior than a single Tikhonov distribution. 

There is a tradeoff in the selection of the threshold Try- It 
should be low enough so that two components associated with 
different phase trajectories are not clustered but high enough so 
that all the components representing the same phase trajectory 
will be clustered together. 



C. Using CMVM to Cluster Mixture Components 

Algorithm 1. selects a set J of components indices from 
mixture d24l >. These components, represent the same phase 
estimate and are approximated using a single Tikhonov. The 
mixture composed of components from (l24l with indices from 
the set J is: 
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Where A is a normalizing constant. 

Using Theorem ( 12. U and skipping the algebraic details, the 
CMVM operator for (ED, is: 
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Similar results apply also for the backward recursion. The 
computational complexity of implementing a modified bessel 
function in a system is prohibitively expensive. We therefore 
present the following approximation 

log(7 (fc)) w k - \ log(fc) - \ log(27r) (36) 
which holds for k > 2, i.e. reasonably narrow distributions. 
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Thus, the simplified approximated versions of (l35l l and (l34i > 



are 
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We also use the approximation for the modified bessel function 
in the computation of a;. 

D. Computation of P u (c k ) 

Another advantage of using the modified canonical model 
comes to effect when we use the forward-backward scheduling 
of the message passing algorithm to compute the message 
P u (c k ). This message describes the LLR of a code symbol 
and the correct estimation of this measure is crucial for the 
decoding of the LDPC. 

We insert our modified canonical model ([19), (l20l i into ([T8~t : 
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Thus P u {c k ) is a weighted summation of four components 
which can be interpreted as conditioning on the probability 
that a phase slip has occurred for each recursion (forward and 
backward). This will ensure that the computation of P u {c k ) is 
based on the most reliable phase posterior estimations, even 
if a phase slip has occurred in a single recursion (forward or 
backward). 

V. Complexity 

As stated in the introduction, DP suffers from a complexity 
limitation. The DP algorithm quantizes the phase symbols and 
performs the SPA on high resolution values of the phase. In 
table 1. we compare the computational load for DP, BARB 
and the algorithm proposed in this paper. 

TABLE I 

Computational load per code symbol per iteration for M-PSK 

CONSTELLATION 
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Fig. 2. Frame error rate - BPSK and wiener phase noise 



We can see that we can achieve a high level of accuracy 
while maintaining a low computational load. Therefore, the 
algorithm proposed in this paper provides an improved tradeoff 
between accuracy and complexity than BARB. 

VI. Numerical Results 

Monte Carlo simulation results for the various algorithms 
are shown in Fig|2] A length 4608 LDPC code with rate 0.889 
was used. The bits were mapped to a BPSK constellation 
and the phase noise model used was a wiener process with 
(Ta = 0.1[rads/symbol]. A single pilot was inserted every 80 
symbols and the parameter To for the selection algorithm was 
chosen to be 2.2. The DP algorithm was simulated using 8 
quantization levels. 

As shown in Fig [2] BARB demonstrates a high error floor. 
This is because of the large phase noise variance and large 
spacing between pilots which causes the posterior estimation 
to become uniform and thus does not provide information 
for the LDPC decoder. The high code rate amplifies this 
problem. The new algorithm has a negligible with respect to 
DP algorithm. 

References 

[1] W. Stark A. Worthen. Unified design of iterative receivers using factor 
graphs. IEEE Transactions on Information Theory, 47:843 -849, February 
2001. 

[2] A. Barbieri G. Colavolpe and G. Caire. Algorithms for iterative decoding 
in the presence of strong phase noise. IEEE Journal on Selected Areas 
in Communications, 23:1748 -1757, September 2005. 

[3] S. Kullback and R. A. Leibler. On information and sufficiency. The 
Annals of Mathematical Statistics, 22:79-86, March 1951. 

[4] KV. Mardia and Jupp P. Directional Statistics. John Wiley and Sons Ltd., 
2000. 



M is the constellation order, L is the number of quantization 
levels and Q is a parameter for the DP algorithm explained in 

na. 



